Ionization front-driven turbulence in the clumpy interstellar medium 
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We present 3D radiation-gasdynamical simulations of an ionization front running into a dense 
clump. In our setup, a BO star irradiates an overdensity which is at a distance of 10 pc and modelled 
as a supercritical 100 Mq Bonnor-Ebert sphere. The radiation from the star heats up the gas and 
creates a shock front that expands into the interstellar medium. The shock compresses the clump 
material while the ionizing radiation heats it up. The outcome of this "cloud-crushing" process is a 
fully turbulent gas in the wake of the clump. In the end, the clump entirely dissolves. We propose 
that this mechanism is very efficient in creating short-living supersonic turbulence in the vicinity of 
massive stars. 

PACS numbers: 47.40. Nm,97.10.Bt,98.38.Am 



I. INTRODUCTION 

Massive stars strongly influence the environment in which they have formed by stellar winds, ionizing radiation and 
supernova explosions. An ionization front which expands into the ambient interstellar medium (ISM) and hits a dense 
clump may compress it so heavily that gravitational collapse might be triggered. On the other hand, ionization heats 
up the material and can photoevaporate the clump. The remaining material of this competition may form low-mass 
stars [lj or brown dwarfs [![. 

The interaction of a shock with a dense clump is called "cloud-crushing". The cloud-crushing scenario has been 
studied numerically both for supernova shocks and ionization fronts. The first extensive studies of the fate of the 
shocked cloud already showed that strong vortex rings can be produced [3]. The mixing properties of the cloud 
depend sensitively on the initial density distribution [4]. Furthermore, simulations of dense clumps exposed to an 
ionizing flux but without strong shocks show the generation of kinetic energy [5] and fragmentation of the clump 
Radiation-gasdynamical simulations have also been used to match observations in H II regions, especially the Eagle 
Nebula @,[|. 

Although ionizing radiation injects a significant amount of energy into the ISM, it does not seem to be an important 
driving mechanism of interstellar turbulence on a global scale [9]. However, the cloud-crushing process generates a 
considerable amount of turbulence locally in the wake of the cloud. We find that the motion of the cloud material is 
mostly supersonic while the ambient gas behind the front moves only subsonically. The continuous heating limits the 
lifetime of the dense material, but the supersonic motions are maintained until the cloud disperses. This is contrary 
to the situation in jet-clump interactions, where the situation is less clear with some studies showing mostly subsonic 
motions [To| while others claim supersonic velocity fields [ITj ]. 

II. NUMERICAL METHODS 

We perform 3D radiation-gasdynamical adaptive mesh simulations with the FLASH code [H], which is based on 
the PARAMESH library [13j]. In addition to the refinement with respect to shocks we also make sure that we resolve 
the Jeans length with a sufficient number of cells to avoid artificial fragmentation [3, EH]. The Jeans length is the 
critical scale for self-gravitating objects, and a dynamical increase of resolution would indicate a runaway collapse. 
However, the run discussed in this paper show collapse only temporarily. 

Additionally to the standard hydrodynamics and self-gravity, we include the ionization feedback from a massive 
star. The radiation feedback is included via a raytracing approach [16j]. There are no radiation pressure terms in the 
Euler equations; in this module coupling between hydrodynamics and radiation takes place only through thermody- 
namics. We solve a rate equation for the ionization fraction including collisional ionization and photoionization as 
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well as radiative recombination. The energy equation contains a photoionization heating rate; the equation of state 
is isothermal. 



III. SETUP 

The computational domain has dimensions (6.0 x 2.4 x 2.4) • 10 19 cm with an effective resolution of 512 x 256 x 256 
cells. We place a self-gravitating supercritical gas sphere with Bonnor-Ebert density profile [13, EH, with mass 
M = 100 Mq in the centre of the domain. It slowly rotates with an angular velocity of ft = 7.83 • 10 -15 rad/s around 
the z-axis corresponding to a ratio of rotational to gravitational energy of (3 = 0.02. Additionally to the rotational 
velocity, a turbulent velocity field with a magnitude of at most 50% of the sound speed is added. The temperature of 
the Bonnor-Ebert sphere is T = 20 K, while the ambient gas has T = 90 K. 

The ionizing source is located at the left hand side of the computational domain at (0.0, 1.2, 1.2) • 10 19 cm. It has a 
temperature of T s = 29, 200 K and a luminosity of L s = 24,000^0, representing a B0 star. The radiation heats the 
interstellar gas to T w 6 • 10 4 K. 

IV. RESULTS 

When the simulation starts, the gas next to the source becomes ionized and heated. An ionization front accompanied 
by a shock expands into the ISM. When the shock hits the clump, the clump material is swept away and thereby 
compressed heavily. The shock leaves behind a fully turbulent gas. The hot ionized gas is being mixed with the cold 
clump material while the radiation heats it up continuously. After the shock front passes the clump, the former clump 
disperses completely. 

The time sequence of this process is as following. The shock touches the outer edge of the clump at t — 0.40 Myr. 
At t = 0.53 Myr, it reaches the clump centre at x = 3 • 10 19 cm, and at t = 0.76 Myr, the front has totally enclosed 
the clump. Of course, the dense material delays the shock, so that the wings of the front can propagate faster. They 
enter the shadow zone and meet in the middle at t = 0.77 Myr. The clash of these wings is an important driving 
mechanism of the turbulence seen in these simulations. It builds up an extended turbulent wake while the ionization 
front propagates further into the ISM. At t = 1.09 Myr, the shock reaches the boundary of the computational domain 
at x = 6 • 10 19 cm. The simulation stops at t = 1.49 Myr when the cloud has dissolved. Some of the stages of this 
time evolution are depicted in figure [H 

V. DISCUSSION 

A. Energy balance 

We start our analysis with a brief look at the energy balance in the flow (see figure [2]) . The plot shows the total 
energy E, internal energy Ei nt and kinetic energy i^kin in erg as a function of time t. The B star provides energy to the 
system by ionising material. This contribution is predominantly transferred into internal energy by the photoionization 
heating, only a small fraction is converted into kinetic energy. For example, at t = 0.5 Myr, the ratio of internal over 
kinetic energy is E[ nt /E^[ n w 57. 

Since the luminosity of the star is constant in time, the energy transferred from the star to the gas in the com- 
putational domain grows linearly. This explains qualitatively the form of E(t) in figure O A quantitative analysis 
is difficult however, since geometric effects have to be taken into account appropriately. One would have to account 
for the fact that the star emits its radiation isotropically, while it is not at the centre of a spherically symmetric 
computational domain, but on one face of a Cartesian box. 

B. Mach numbers 

The simulation shows that the cloud-crushing flow is largely dominated by supersonic motion. This is because the 
cloud material is cold, so that the sound speed is much lower than in the hot gas behind the ionization front, where 
a wind with M ~ 0.2-0.4 is observed. The cause of the wind is that the photoionization heating is stronger close the 
source, which leads to a pressure gradient and a corresponding flow. Since the wind prevails in the largest part of the 
domain, namely the hot post shock gas, the mean Mach number .M m ean is always below unity, while the maximum 
Mach number M. max , which is reached at crushing, can be greater than 6. 
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Figure [3] depicts 7W m ax and A4 mean as a function of time. The maximum Mach number traces the shock ahead of 
the ionization front. Within a time of 0.15 Myr, the shock accelerates up to a constant velocity of M = 4. For a 
short time of another 0.15Myr the velocity seems to saturate, but the continuous photoionization heating accelerates 
the shock again up to M = 6.5. At this point, when M m &x is maximal, the shock collides with the dense clump. As 
it hits the high-density gas, the shock front moves more slowly. The gas decelerates, so that A^max decreases again. 
However, the wings of the shock that were not affected by the clump can enter the shadow zone, which leads to a peak 
in Almax after 0.7 Myr. Then these wings collide, which stops the motion in ^/-direction, so that the Mach number 
decreases further. But the collision of the wings also leads to an acceleration in positive x-direction, which can be 
seen in a series of peaks from 0.8 Myr to 1.1 Myr. After 1.1 Myr, the shock leaves the computational domain, resulting 
in a sharp drop in M mdiX . This demonstrates that the highest Mach numbers are only reached in the shock front 
itself, not in the shock-generated turbulence behind the front. The motion in the turbulent wake is mostly supersonic 
with Ai below 3. The mean Mach number grows continuously until the shock leaves the domain, whereafter .Mmean 
declines slowly. The heating of the gas does not change .Mmean. 

The different stages of the simulation can also be recognized in the probability density functions (PDFs) of the 
Mach number M. Figure [4] shows mass- weighted PDFs at the moment of cloud-crushing (t = 0.49 Myr) and after the 
shock has left the domain (t = 1.16 Myr to t = 1.45 Myr). At the crushing time, the high M above 2 all belong to 
the shock. The shock then excites supersonic turbulence in the wake, but away from the shock M above 3 is very 
rare. While most of the dense gas is supersonic, most of the domain is dominated by low Mach number flows, both 
at crushing time and afterwards. 

The PDFs of the total Mach number M should be compared with the Mach number given only by the turbulent 

velocity fluctuations, v t urb = yj^y + ^ ? which is ^turb/c s , where c s is the local speed of sound. These plots are shown 

in figure [H Since the bulk motion in ^-direction is no longer taken into account, the Mach numbers are significantly 
lower. At the moment of cloud-crushing, the turbulent Mach number is only slightly supersonic, while it is totally 
subsonic afterwards. Hence, the bulk motion of the shock (and also the transport of momentum by the wind) is 
important to reach the high Mach numbers observed above. 

We are also interested in the fraction of mass which moves supersonically. In figure [6] we plot the ratio of the mass 
of supersonic gas M sup and the total gas mass in the computational domain M tot . Despite of the complicated mixing 
processes, M sup /M tot grows roughly linearly until the shock leaves the domain. It is surprising that a significant part 
of the gas moves supersonically, altough most of the energy input is converted into internal energy (see figure [2]) . 
When the shock front leaves the domain, more than 60% of the gas in our domain is in supersonic motion. 

C. Clump structure after crushing 

In figure El we enlarge a part of figure []] belonging to the snapshot at t = 1.27 Myr and show additionally to the 
mass density also the temperature and the velocity components v x and v y . The remains of the dense core have a 
temperature around T w 300 K, while the ambient ionized medium is at T w 3.5 • 10 4 K. The high temperatures in 
the environment give rise to the "rocket effect", which accelerates the gas in positive x-direction [13]. This is because 
the cold gas at the surface of the clump facing the star becomes heated. Thus, it expands into the postshock medium, 
carrying momentum with it, and consequently the clump accelerates. 

The cloud-crushing leads to vortical structures in the wake as can be seen from the lower plots in figure The two 
largest structures around x w 4.2 • 10 19 cm and y w 0.7 • 10 19 cm or y w 1.7 • 10 19 cm, respectively, even have slightly 
negative v x and so does a region around the tip of the former core at x = 3.3 • 10 19 cm and y = 1.2 • 10 19 cm. Averaged 
over the whole volume, however, the wind, wich moves with v x ~ 1.5 • 10 6 cm/s, causes a bulk motion of the gas in 
positive x-direction. In order to measure the turbulent components of the velocity field, it is better to focus on the 
transversal directions. The peak amplitude of the velocity component in ^/-direction, for example, is about half of the 
maximum velocity in ^-direction. The large vortical structures discussed above have total velocities around 10 6 cm/s. 
The upper vortex rotates clockwise, while the lower vortex rotates counter-clockwise. 

D. Line-of-sight velocity profiles 

A connection to observations of molecular clouds can be made by looking at velocity profiles along a certain line- 
of-sight [20]. Examplary, we show in figure [8] a profile for v y of the dense core after crushing at time t = 1.27 Myr 
(compare also figure El). The width of the beam is 3.0- 10 19 cm < x < 4.0- 10 19 cm and 0.7- 10 19 cm < z < 1.7- 10 19 cm, 
while y ranges over the whole box width. Note that the profile is a mass-weighted histogram, not a normalized PDF. 
We do the same calculation again, but this time only considering gas that has T < 10 3 K (low T case) or T > 10 4 K 
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FIG. 1: The main stages of the simulation are shown as 2D cuts of the logarithmic mass density log(p / gem -3 ) in the midplane. 
The different snapshots are taken at times t = 0.16 Myr, t = 0.55 Myr and t = 1.27 Myr, respectively. The turbulent wake 
behind the former clump at the last stage is clearly visible. 



(high T case). From the figure we see that the high velocity tails are exclusively related to the high temperature gas. 
The peak at low velocities comes mainly from both the low temperature medium as well as a warm envelope with 
temperatures between the cuts 10 3 K <T< 10 4 K. Additionally to the histograms, we have fitted Gaussians to the 
low and high T data. Their variance gives a measure for the turbulent velocity and can be related to the width of 
spectral lines that are being oberserved. 



VI. CONCLUSION 



We have seen that cloud-crushing by ionization fronts can lead to short-living supersonic turbulence. Altough only 
a minute fraction of the energy input is converted into kinetic energy, up to 60% of the affected gas is supersonic. 
While it is mainly the cold gas that is highly supersonic, it is the hot gas that moves the fastest. The bulk motion of 
the shock is an important contribution to the supersonic flow, since the transversal fluctuations are at best slightly 
supersonic. 
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FIG. 2: The total energy E, internal energy Ei n t and kinetic energy i^kin are shown as function of time t. The major contribution 
to E comes evidently from E- in t, while i^kin is between one and two orders of magnitude smaller. The point in time where the 
shock, which carries most of the kinetic energy, leaves the computational domain can be distinguished easily. 
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FIG. 3: The average and maximum Mach numbers in the flow provide information on the dominance of shocks in the flow. To 
compare these two numbers, the mean Mach number A^mean is displayed amplified by a factor of 5. The peaks in .Mmax can 
be associated with events in the cloud-crushing scenario. 
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FIG. 4: The mass-weighted PDFs of the Mach number M between times t — 0.49 Myr and t — 1.45 Myr. The shock front 
shows up as a peak in the regime of high Mach numbers for the first point in time. At t — 1.16 Myr, the shock has left the 
computational domain, and the turbulence with its supersonic Mach numbers decays. 




FIG. 5: The mass- weighted PDFs of the "turbulent Mach number" ^turb/c s with v t urb = y/vy + vl for the same points in time 
as in figure [4] At cloud-crushing, the turbulent field ^turb is slightly supersonic, while the flow in the wake is only subsonic. 
This shows that the bulk motion contributes significantly to the high total Mach numbers. 
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